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We study precise calculation of neutralino relic density in the minimal supergravity 
model. We compare the exact formula for the thermal average of the neutralino 
annihilation cross section times relative velocity with the expansion formula in 
terms of the temperature, including all the contributions to neutralino annihilation 
cross section. We confirm that the expansion formula fails badly near s-channel 
poles. We show that the expansion method causes 5-10 % error even far away from 
the poles. 



1 Introduction 

In recent years, observations of cosmic microwave background radiation have 
been providing experimental values for cosmological parameters with some 
accuracy. Current experimental bound on the fractional energy density of cold 
dark matter is ||,|| 

f^CDM/l' - 0.1 - 0.3, (1) 

where h « 0.7 is the parameter in the Hubble constant Ho = 100 h km/sec/Mpc 
iQ. This bound is expected to be improved in the near future, so precise 
calculation of the relic density becomes important. 

Particle physics should provide a particle-theoretical explanation for dark 
matter. The minimal supergravity model |^ is one of the most promising 
candidate for new physics beyond the standard model. In this model, a discrete 
symmetry, so-called R-parity, is introduced to avoid rapid proton decay. One of 
the consequences of introducing the R-parity is that the lightest superparticle 
(LSP) is stable, and the LSP becomes a good candidate for cold dark matter Q. 
The LSP in this model is mostly the lightest neutralino x which is a mass 
eigenstate given by a linear combination of neutral gauginos and Higgsinos 

X = NiiB + Ni2W^ + Ni3H° + Ni4H2- (2) 

In general, supersymmetric models have a huge number of unknown pa- 
rameters. However, the minimal supergravity model has only a small number 
of parameters, so it has reasonable predictive power. There are five unknown 
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parameters in this model 

mo, mi/2, A, tan/3, sgn(^), (3) 

where mo, mj^/2 ^nd A represent a common scalar mass, a common gaugino 
mass, and a common trilinear scalar coupling, respectively. They parametrize 
soft supersymmetry breaking terms at the GUT scale « 2 x 10^^ GeV. tan (3 
is the ratio of vacuum expectation values of the two neutral Higgs fields. ^ 
denotes the Higgs mixing mass parameter. In the case of the minimal super- 
gravity model, the absolute value of /i is determined from the condition for 
consistent radiative electroweak symmetry breaking, and only the sign is a 
free parameter. In this work, we assume that there are no CP violating phases 
in these parameters (||) for simplicity. 

2 Calculation of the relic density 

In this section, we first briefly review the standard calculation of the neutralino 
relic density in the minimal supergravity model. After that, we explain what 
is necessary to calculate it precisely. 

We evaluate the relic density at present, starting from thermal equilibrium 
in the early universe. The time evolution of the neutralino number density in 
the expanding universe is described by the Boltzmann equation 

^ + 3Hn^ = - (a«M0i> (n^ - n^^^) , (4) 

where H is the Hubble expansion rate, describes the actual number density 
of the neutralino, while n'^ is the number density which the neutralino would 
have in thermal equilibrium, a denotes the cross section of the neutralino pair 
annihilation into ordinary particles. VM01 is so-called M0ller velocity which 
can be identified as the relative velocity between the two neutralinos. (cruM0i) 
represents the thermal average of aVMai- 

In the early universe, the neutralino is in thermal equilibrium = n'^. As 
the universe expands, the neutralino annihilation process freezes out, and after 
that the number of the neutralinos in a comoving volume remains constant. 

In the annihilation cross section a in eq. (|^) , there are a lot of final states 
which contribute: XX — *■ //; hh, WW, ZZ, Zh, etc. Among these final 
states, fermion pairs // usually give dominant contributions. 

Using an approximate solution to the Boltzmann equation (^), the relic 
density = m^n^ at present is given by 

1 fT^\^ ^ 1 
Mpi V / ^ Jq" dx{avMa\) ' 
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where x ~ T/m-^ is a temperature of the neutraUno normahzed by its mass. 
T-^ and Tj are the present temperatures of the ncutralino and the photon, 
respectively. The suppression factor (T^/T^Y ~ 1/20 foUows from the entropy 
conservation in a comoving volume [gj. A/pi denotes the Planck mass, xp ~ 
1/20 is the value of x at freeze-out, and is obtained by solving the following 
equation iteratively: 



where Gat is the Newton's constant, and 5, represents the effective number of 
degrees of freedom at freeze-out. 

The relativistic formula |^ for thermal average in eq. (^) can be written 



as an integration over one variable |10 
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{crvM0i) = o^4^r^2^_ /^^^ / dscr{s){s - ^m\)^sKx ( — ) , (7) 



where Ki (i = 1,2) arc the modified Bessel functions. The cross section a{s) 
is a complicated function of s in general, so we have to evaluate the above 
integration numerically to calculate the thermal average. 

If the normalized temperature x is small enough, we may be able to use 
an expansion formula in terms of x for the thermal average (|^) 

(CTWM01) = a + bx. (8) 

This expansion, neglecting the higher order of x, is widely used in literatures. 
In the case that the neutralino is the LSP, the a coefficient for a fermion pair 
production is proportional to the fermion mass due to the Majorana nature 
of the neutralino. On the other hand the b coefficient for the same final state 
includes a contribution which is not proportional to the fermion mass, so the b 
coefficient is much larger than the corresponding a coefficient for each fermionic 
final state: a < 0. 

In order to calculate the relic density precisely enough, we have to take 
into account the following items: 

(i) all the contributions to the annihilation cross section 

(ii) the exact formula for the thermal average 

(iii) solving the Boltzmann equation exactly 

(iv) coannihilation (This becomes important when the mass of the next LSP 
is close to the LSP mass |l|-|l|.) 
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In this work, however, we concentrate on (i) and (ii) to compare the expan- 
sion formula with the exact one, and we don't take into account (iii) and (iv) 
for simphcity. Namely, we use the approximate solution (||) to the Boltzmann 
equation, and we neglect coannihilation effects. 

We have derived analytic expressions for the exact annihilation cross sec- 
tion and have obtained a and b coefficients analytically for every contribution 
including interference terms. These interference terms are neglected in most 
literatures, but we found that they can give significant contributions in some 
cases. Some of the analytic expressions can be found in literatures for both 
the exact cross section and the expansion coefficients JT8|,[r^. 

3 Thermal average — exact formula vs. expansion 

Let us see the rough behavior of the integrand in eq. @. In the case of interest, 
the argument of the function Ki is much larger than unity, since T ^ m^/20 
and ^/s > 2m^. Therefore the thermal average can be written as a convolution 
of the cross section with a function which decays exponentially: 

/•OO 

(crvMai) ~ / ds(T{s)F{s), (9) 

J 4)n| 

where the function F{s) has an exponential suppression factor as 

F{s) (X e-^/^. (10) 

Naively, it seems that the expansion should converge quickly, since the function 
F{s) in eq. decays exponentially as s increases [||. However this is not 
always true. It is known that the expansion method fails badly when the 
aimihilation cross section changes rapidly with s ||l^, |l^, |l^ . This happens, 
e.g, near s-channel poles and thresholds of new channels. In the following we 
compare the expansion result with the exact one, including all the contributions 
to the annihilation cross section. 

4 Numerical results 

The spectrum and the couplings in the minimal supergravity model at the 
weak scale are calculated by solving renormalization group equations with 
boundary conditions at the GUT scale and implementing radiative electroweak 
symmetry breaking. In this analysis, we use an existing code "SUSPECT" [po| 
to calculate the spectrum and the couplings. 

In the following, we present the results of our calculation. Fig. |l| shows 
the integration of the thermal average J^'' dx (avu^i) as a function of mi/2 
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for tan/3 = 2, too = 200 GeV, A — 0, ii > 0. Note that the rehc density in 
eq. (||) is proportional to the inverse of this integration. The sohd Hne and the 
dashed hne represent the exact resuh and the expansion resuh, respectively. 

We have ploted the total contributions and the individual contributions 
for relevant final states. The exact results show similar behaviors with those 
in a previous analysis 

There are two peaks in the total contributions. One is from the Z boson 
contribution, and it takes place when 2711-^ ~ mz is satisfied. The other is from 
the lightest Higgs (h) contribution, and it happens when 2to;^ « ruh- Around 
the poles, we see huge difference between the exact result and the expansion. 
Away from the poles, the difference is not so huge. 

We plot the relic density vs. toi/2 in Fig- H (e^) for the same parameters 
as those in Fig. |^. In the region of the poles, we see a dip of the relic density. 
Away from the poles, the relic density is too large to satisfy the experimental 
constraint (pj). 

In Fig. I (b), we plot the ratio Jlcxpansion/^^ exact; where ^c-^act (^expansion) 
denotes the relic density calculated with the exact (expansion) formula for the 
thermal average. In the pole regions, the expansion is quite different from the 
exact one by a large factor . Also we find a difference of about 5 % in this 
ratio even far away from the pole. 

In a very large tan /3 case, the result has a different feature. In Fig. |^ (a) 
and (b), we show similar results for tan/3 = 50, too = 600 GeV, yl = 0, /i > 
0. The relic density is small due to enhancements of the pseudoscalar Higgs 
exchange contribution and the heavier neutral Higgs exchange contribution in 
the bb final state for a large tan /3. The are two reasons for these enhancements. 
First, couplings of the pseudoscalar Higgs and the heavier neutral Higgs to the 
bottom quark become large for a large tan /3. Secondly, masses of heavy Higgses 
become smaller in a large tan/3 region. Because of these enhancements in the 
annihilation cross section, the relic density in Fig. ^ (a) becomes small enough 
to satisfy the experimental constraint (|l|) for a wide region of 600 GeV < mi/2 
< 800 GeV. As for the ratio of the expansion result to the exact one in Fig. || 
(b), we see a difference more than 10 % for 600 GcV < toi/2 < 800 GeV. Note 
that we find this relatively large difference in the interesting region where the 
relic density satisfies the experimental constraint (|l|). 

We show a contour plot for the quantity 2to^ — rriA for tan/3 — 50, A = 
and fj, > in Fig. ^ where ruA denotes the pseudoscalar Higgs mass. Along 
the thin solid line, this quantity vanishes so that there are enhancements in the 
annihilation cross section from the pseudoscalar Higgs pole contribution and 
the heavier neutral Higgs pole contribution through the bb final state. Fig. ^ 
shows that the region 600 GeV < TOi/2 < 800 GeV in Fig. | (b), where the 
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deviation of ficxpansion from Qcxact is about 10 %, differs from the pole region 
(i.e. the zero of the quantity 2to^ — ttia) in mo by about 150 GeV. Despite of 
this difference, we still have a sizable error of 10 % in the expansion result, as 
we have seen in Fig. || (b) . 

In Fig. H and Fig. ^, we find 5-10 % difference even away from the pole. 
Thus the expansion never succeeds to approximate the exact result within the 
accuracy of 5 %. This implies that the coefficient of the next order cx^ in the 
expansion (H) has the same order with the b coefficient c/b w 0(1). 

5 Conclusions 

We have presented results of precise calculation of neutralino relic density in 
the minimal supergravity model. We have calculated all the contributions 
to neutralino annihilation cross section analytically. We have compared the 
exact formula for the thermal average of the annihilation cross section times 
relative velocity with the expansion formula in terms of the temperature. We 
have confirmed that the expansion formula fails badly near s-channel poles. 
We have shown that the expansion method causes 5-10 % error even far away 
from the poles. 
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Figure 1: Integration of the thermal average J^'' dx (crvyi^i) as a function of rrii/^ fo^' tan /3 
= 2, mo = 200 GeV, A = 0, ij, > 0. 
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(a) (b) 

Figure 2: The relic density flh^ (a) and the ratio Ocxpansion/f^oxact (b) for tan/3 = 2, mo 
= 200 GeV, A = 0, M > 0. 




(a) (b) 

Figure 3: Similar results to Fig. for tan/3 = 50, mo = 600 GeV, A = 0, /t > 0. 
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tan/S = 50, A = 0, fj, > 
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Figure 4: A contour plot of 2m-^ — niA in the mo-mx/2 plane for tan /3 = 50, A = and fi > 
0. The region 'Charged LSP' is excluded because the lighter stau is the LSP in this region. 
In the region 'No EWSB', the electroweak symmetry breaJjing does not occur. 
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